#include<iostream>
using namespace std;
#include<cmath>

/*
 * Naive method to calculate the derivative of func at x,
 * using formula:
 * 	f'(x) = (f(x+h) - f(x-h))/2h + O(h^2)
 * h is the step, note that if it's too small, the round off error becomes large enought that the derivative makes no sense.
 */
double naive_deriv( double x, double h, double (*func)(double) ){

	if( fabs(h)<1E-9 ){	
		cout<<" Error in deriv_f(): h <= 1E-9 "<<endl;
		exit(1);
	}
	return ( func(x+h) - func(x-h) )/(2*h);
}

/*
 * Richardson算法：
 * 首先计算
 * 	D^k_0 = (f(x+h/2^k) - f(x-h/2^k))/ ( h/2^{k-1} ), k=0,1,...,n
 * 然后根据递推公式
 * 	D^k_m = D^{k+1}_{m-1} + \frac{ D^{k+1}_{m-1} - D^k_{m-1} }{ 4^m -1 }
 * 逐步推得　D^0_n 作为返回值，其误差为　O(h^{2(n+1)}
 *
 * 例：
 *
 * D^0_0
 * D^1_0	D^0_1
 * D^2_0	D^1_1	D^0_2
 *
 * 我们定义一个　3x3　的矩阵 deriv，然后逐列计算矩阵元，最后返回　D^0_2, 即 deriv[2][2] 的值。
 * 
 * n 不能是负数，否则程序终止
 */
double Richardson_double_deriv(double x, double h, int n, double (*func)(double) ){

	double coef, y;

	double ** deriv = new double * [ n+1 ];//开辟动态内存，得到新的 (n+1)*(n+1) 矩阵
	for(int i=0;i<n+1;i++) deriv[i] = new double [n+1];
	
	for(int i=0;i<=n;i++){ // 赋值第一列 D^k_0, k=0,1,...,n
		coef = pow(0.5,i);// 1/(2^i)
		deriv[i][0] = ( func(x + h*coef ) + func(x - h*coef) - 2* func(x) ) / (h* coef) / (h*coef);
	}
	for(int i=1;i<=n;i++){ // 赋值第　1,...,n 列
		for(int j=i;j<=n;j++){ // j 是行数，i是列数
			deriv[j][i] = deriv[j][i-1] + (deriv[j][i-1] - deriv[j-1][i-1])/ (pow(4.0,i) -1 );//　这一行与递推公式对应
		}
	}
	//cout<<" error is O(h^"<<2*n+2<<")\n";
	y = deriv[n][n];
	for(int i=0;i<=n;i++)delete [] deriv[i];//释放动态内存
	delete [] deriv;

	return y;
}

int main(){

	double x=10;
	double h = 0.1;
	printf("f''(10) = %.15lf, error = %.15lf\n",
		Richardson_double_deriv(x, 0.1, 3, exp),
		- exp(x) + Richardson_double_deriv(x, 0.1, 3, exp) );
	return 0;
}
